﻿using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using OxyPlot;
using WaveletStudio;
using WaveletStudio.FFT;
using WaveletStudio.Functions;
using WaveletStudio.Wavelet;

namespace AutoCorrelatorGUI.Model
{
   public class Well_Transform
   {
      public Well_Transform()
      {
         FourierData = new List<DataPoint>();
         WaveletData = new List<DataPoint>();
         PowerData = new List<DataPoint>();
      }

      public List<DataPoint> FourierData { get; private set; }
      public List<DataPoint> WaveletData { get; private set; }
      public List<DataPoint> PowerData{get; private set;}

      private uint NextPowerOf2(uint n)
      {
         n--;
         n |= n >> 1;   // Divide by 2^k for consecutive doublings of k up to 32,
         n |= n >> 2;   // and then or the results.
         n |= n >> 4;
         n |= n >> 8;
         n |= n >> 16;
         n++;

         return n;
      }

      private MotherWavelet GetMotherFromEnum(MotherWaveletTypes type)
      {
         switch (type)
         {
         case MotherWaveletTypes.Coiflet1:
            return CommonMotherWavelets.GetWaveletFromName("coif1");
         case MotherWaveletTypes.Coiflet2:
            return CommonMotherWavelets.GetWaveletFromName("coif2");
         case MotherWaveletTypes.Coiflet3:
            return CommonMotherWavelets.GetWaveletFromName("coif3");
         case MotherWaveletTypes.Coiflet4:
            return CommonMotherWavelets.GetWaveletFromName("coif4");
         case MotherWaveletTypes.Coiflet5:
            return CommonMotherWavelets.GetWaveletFromName("coif5");
         case MotherWaveletTypes.Daubechies2:
            return CommonMotherWavelets.GetWaveletFromName("db2");
         case MotherWaveletTypes.Daubechies3:
            return CommonMotherWavelets.GetWaveletFromName("db3");
         case MotherWaveletTypes.Daubechies4:
            return CommonMotherWavelets.GetWaveletFromName("db4");
         case MotherWaveletTypes.Daubechies5:
            return CommonMotherWavelets.GetWaveletFromName("db5");
         case MotherWaveletTypes.Daubechies6:
            return CommonMotherWavelets.GetWaveletFromName("db6");
         case MotherWaveletTypes.Daubechies7:
            return CommonMotherWavelets.GetWaveletFromName("db7");
         case MotherWaveletTypes.Daubechies8:
            return CommonMotherWavelets.GetWaveletFromName("db8");
         case MotherWaveletTypes.Daubechies9:
            return CommonMotherWavelets.GetWaveletFromName("db9");
         case MotherWaveletTypes.Daubechies10:
            return CommonMotherWavelets.GetWaveletFromName("db10");
         case MotherWaveletTypes.Haar:
            return CommonMotherWavelets.GetWaveletFromName("haar");
         case MotherWaveletTypes.DiscreeteMeyer:
            return CommonMotherWavelets.GetWaveletFromName("dmeyer");
         case MotherWaveletTypes.Symlet2:
            return CommonMotherWavelets.GetWaveletFromName("sym2");
         case MotherWaveletTypes.Symlet3:
            return CommonMotherWavelets.GetWaveletFromName("sym3");
         case MotherWaveletTypes.Symlet4:
            return CommonMotherWavelets.GetWaveletFromName("sym4");
         case MotherWaveletTypes.Symlet5:
            return CommonMotherWavelets.GetWaveletFromName("sym5");
         case MotherWaveletTypes.Symlet6:
            return CommonMotherWavelets.GetWaveletFromName("sym6");
         case MotherWaveletTypes.Symlet7:
            return CommonMotherWavelets.GetWaveletFromName("sym7");
         case MotherWaveletTypes.Symlet8:
            return CommonMotherWavelets.GetWaveletFromName("sym8");
         default:
            Debug.Fail("Unknown Mother Wavelet Type");
            return null;
         }
      }

      public void Calc(Well_Output output, double intervalRateMS, MotherWaveletTypes motherWavelet, uint maxDataCount)
      {
         uint N = (uint)output.Data.Count;
         uint next2 = NextPowerOf2(maxDataCount);

         int pow2Count = (int)Math.Log(next2, 2.0);

         uint zeroPad = next2 - N;

         //build the data necessary for the fft padding zero on either end
         double[] data = new double[next2];

         for (int i = 0; i < N; i++)
         {
            data[i + zeroPad / 2] = output.Data[i].Y;
         }

         Signal dwtSignal = new Signal(data, (int)Math.Round(1.0 / (intervalRateMS / 1000.0)));
         dwtSignal.Finish = (intervalRateMS / 1000.0) * (next2 - 1);

         var dwtOutput = DWT.ExecuteDWT(dwtSignal, GetMotherFromEnum(motherWavelet), pow2Count, SignalExtension.ExtensionMode.ZeroPadding);

         List<double> finalDWT = new List<double>();

         for (int i = dwtOutput.Count - 1; i >= 0; i--)
         {
            finalDWT.AddRange(dwtOutput[i].Details);
         }

         //Get the max value
         double maxDWT = 0.0;

         for (int i = 0; i < finalDWT.Count; i++)
         {
            maxDWT = Math.Max(Math.Abs(finalDWT[i]), maxDWT);
         }

         WaveletData.Clear();

         for (int i = 0; i < finalDWT.Count; i++)
         {
            WaveletData.Add(new DataPoint(i, finalDWT[i] / maxDWT));
         }

         //Puts a 0 for every odd index which is the complex component
         data = WaveMath.UpSample(data, 2, false);

         ManagedFFT.Instance.RealFFT(ref data, true);

         FourierData.Clear();
         PowerData.Clear();

         double samplingRate = 1.0 / (intervalRateMS / 1000.0);

         double maxMag = 0.0;

         double[] outMags = new double[next2 / 2];

         for (int i = 1; i <= (next2 / 2); i++)
         {
            double mag = Math.Sqrt(data[2 * i] * data[2 * i] + data[2 * i + 1] * data[2 * i + 1]);

            maxMag = Math.Max(maxMag, mag);

            outMags[i - 1] = mag;
         }

         for (int i = 0; i < outMags.Length; i++)
         {
            double f = (i * samplingRate / 2.0) / (next2 / 2.0); //Find the middle point for the frequency

            double dB = 20.0 * Math.Log10(outMags[i] / maxMag);

            double power = (1.0 / (double)next2) * outMags[i] * outMags[i];

            FourierData.Add(new DataPoint(f, dB));
            PowerData.Add(new DataPoint(f, power));
         }
      }
   }
}
